IMPLEMENTING UNSUPERVISED LEARNING FOR SPLITTING COUNTRIES IN THE WORLD IN TERMS OF POSITION AGAINST SPREAD EFFECT

1. INTRODUCTION

Este trabajo es una continuación de este, por lo tanto lo primero que tendremos que tener claro, para poder interpretar los diferentes clusters que vayamos sacando y tener un criterio aparte de los matemáticos para quedarnos con uno u otro tipo de cluster, será recordar cuál era el objetivo para el que se recogieron las variables que se utilizaron en esa primera parte del proyecto y que se utilizarán también en esta. El objetivo era, en definitiva, ver cómo se posicionaban los países en función del spread y el bounty.

Históricamente se ha hablado siempre del primer, segundo y tercer mundo; sin embargo, en los últimos años se habla más de países desarrollados (bien posicionados respecto al efecto de spread generado por el progreso tecnológico) y países en vías de desarrollo (los que se están quedando atrás con respecto al resto debido a este progreso tecnológico, ya que se están comiendo una parte más pequeña de la tarta que este genera).

Por lo tanto, estaremos cómodos tanto con 2 como con 3 clusters. Es posible que encontremos alguna composición de un número mayor de clusters con la que nos sintamos cómodos porque al analizarla tenga bastante sentido y una interpretabilidad suficiente. Sin embargo, a priori buscaremos encontrar 2 o 3 clusters.

Como todos los gráficos de cluster los vamos a realizar sobre los dos primeros componentes principales, remito al lector a este link, en el que se especifica y se analizan los primeros componentes principales.

Adicionalmente, incluimos un pequeño párrafo con las variables que son objeto de estudio:

Gdp Per Capita: la sacamos de WB
Inflation: la sacamos de WB. Medida como índice de precios al consumo.
Gni Per Capita: la sacamos de WB
Urban Population: porcentaje de población urbana. La sacamos de UN.
Employment Agriculture: porcentaje de empleo en agricultura. UN
Economy Agriculture: Importancia de la agricultura en la economía. UN.
Urban Population Growth: de 2016 a 2017. UN.
Credit Information: Del 1 al 8 la información de crédito que tiene la gente de un país. WB.
Productivity: GDP/Hour worked. WB.
Renewable Energy: Uso de energía renovable. WB.
Renewable Electricity: Porcentaje de la energía eléctrica que es renovable. WB.
Electricity 2016: WB.
Clean Fuels 2016: Acceso a “clean fuels”. WB.
Infant Mortality Rate: WHO.
Sanitation Services: WHO.
Deaths for Household Pollution: WHO.
Water: WHO.
Fertility Rate: UN.
Account Ownership: WB.
Life Expectancy: UN.
Vulnerable Employment: WB.

La estructura de este report será la siguiente:

Introducción.
Visualización Multivariante.
K-Means Clustering.
PAM.
Model-Based Clustering.
Hierarchical Clustering.
Comparación de Resultados.
Conclusiones.

Eliminamos la primera columna, ya que sólo contiene el índice. También eliminamos la observación en la que la inflación era máxima, pues, como se vio en el tratamiento previo de las variables, que se puede encontrar aqui, esta observación era un outlier.

data$X <- NULL
data <- data[!data$inflation == max(data$inflation), ]

En este caso no realizamos un análisis exploratorio de las variables, porque eso ya lo hicimos en el link que aparece arriba; por lo tanto iremos directamente a probar diferentes métodos de cluster, y buscaremos el que más nos convenza a la hora de separar los países en estos diferentes grupos.

Nos quedamos sólo con las columnas numéricas para la matriz X.

countries <- data[, "country"]

X_unscaled <- data[, -c(22, 23, 24)]


X <- X_unscaled

2. Visualización Multivariante.

Antes de entrar a implementar las técnicas de aprendizaje no supervisado con las que trataremos de agrupar los países, vamos a visualizar de forma sencilla la matriz de correlaciones de las variables, así como los histogramas de la distribución de estas variables. Esto último es lo primero que vamos a hacer, pues así podremos ver qué transformaciones en las variables podríamos hacer para tratar de asemejarlas más a una normal. En este link vemos formas típicas de conseguir que distribuciones asimétricas tomen una forma más parecida a la normal. Además de las que vienen en el enlace, probaremos otras para quedarnos con las que más nos guste, siempre usando el histograma de la variable transformada como criterio. Sine embargo, con el fin de no ensuciar el report, mostraremos directamente sólo el código de transformación de estas variables, sin dibujar de nuevo estos histogramas.

for (i in 1:ncol(data.frame(X))) {
  
  hist(data.frame(X)[ , i], col = "blue")
  
}

En este chunk se pueden ver las transformaciones que les hacemos a las variables para tratar de que sus distribuciones se parezcan lo máximo posible a la Normal. Esto será importante especialmente a la hora de interpretar las probabilidades de los Mixture Models, que probaremos más adelante.

Transformaciones.

X[ , 1] <- X[ , 1]**(1/3)

#X[ , 2] <- X[ , 2]**(1/3)

X[ , 3] <- X[ , 3]**(1/4)

X[ , 5] <- X[ , 5]**(1/3)

X[ , 8] <- X[ , 8]**(1/2)

X[,9] <- X[ , 9]**(2)**(1/3)

X[, 11] <- X[,11]**(1/3)

X[ , 13] <- X[,13]**(1/2)

X[, 15] <- X[,15]**(1/4)

X[,16] <- X[,16]**(1/2)

#X[,17] <- X[,17]**(1/4)

X[,19] <- X[,19]**(1/3)


X <- data.frame(scale(X))

library(GGally)

ggpairs(X)

cor_matrix <- cor(X)

ggcorr(X, label = T, label_round = 2, label_alpha = T)

Podemos ver que la mayoría de las variables tienen bastante correlación, tanto positiva como negativa, con el resto. Esto es buena señal, pues indica que, pares a pares (es decir, en dimensión 2), nuestras variables están bastante relacionadas; esto hará que los modelos de aprendizaje no supervisado que implementaremos a continuación tengan mayor facilidad para separar a los países en grupos.

3. K-MEANS CLUSTERING

Using Euclidean Distance

fit = kmeans(X, centers=5, nstart=100)
groups = fit$cluster
#groups
barplot(table(groups), col="blue")

centers=fit$centers
centers

##   gdp_per_cap  inflation gni_per_cap    urban_p productivity  credit_inf
## 1  0.02419864  0.3327386  0.08761484  0.2245190    0.3409200  0.01201352
## 2 -0.74670445 -0.1864966 -0.79203840 -0.8211653   -0.8148517 -0.61891526
## 3 -1.11412294  0.9075402 -1.18388174 -0.9423741   -1.3117730 -0.94862434
## 4 -0.32684608 -0.2665991 -0.28436011 -0.1384272   -0.2125053  0.47140805
## 5  1.40617124 -0.5015226  1.38975542  1.0172553    1.2224013  0.54036621
##      renew_en renew_electr electricity_2016 clean_fuels_2016 inf_mort_rate
## 1 -0.87348254  -0.78239308        0.5818854        0.7360304    -0.1981290
## 2  0.30659661   0.10239740       -0.6629780       -1.0263136     0.9862122
## 3  1.59717242   0.32985108       -1.8265185       -1.5472748     1.3701156
## 4 -0.02238447   0.44398477        0.5327119        0.1667835    -0.1420763
## 5 -0.48843966  -0.04074668        0.6239364        0.8522864    -1.1611765
##   sanitation_services deaths_household_pollution      water  fert_rate
## 1           0.5867365                 -0.5722864  0.4844916 -0.2879624
## 2          -0.9095397                  0.9925236 -0.6274547  0.7596367
## 3          -1.6589471                  1.1896931 -1.7622550  1.5181949
## 4           0.3162657                  0.2922737  0.3766224 -0.2978457
## 5           0.8436136                 -1.1339905  0.7753450 -0.9189310
##     life_exp urban_pop_growth      emp_agr    eco_agr vulner_empl
## 1  0.2473097       -0.3302692 -0.539902591 -0.4052800 -0.52012687
## 2 -0.8223502        0.5529592  0.700815723  0.6052136  0.69963044
## 3 -1.5232020        1.2979051  1.517756967  1.3075144  1.54357537
## 4  0.2542822       -0.2868240  0.002811144  0.2569874  0.04609847
## 5  1.0281593       -0.6359001 -0.948799884 -1.0873899 -1.01792821
##       acc_own
## 1  0.06864725
## 2 -0.75933271
## 3 -0.91305480
## 4 -0.28961642
## 5  1.21761600

# clusplot
fviz_cluster(fit, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

# Silhouette plot
# The silhouette value in [-1,1] measures the similarity (cohesion) of a data point to its cluster relative to other clusters (separation). 
# Silhouette plots rely on a distance metric and suggest that the data matches its own cluster well.
# The larger the silhouette widths, the better.
d <- dist(X, method="euclidean")  
sil = silhouette(groups, d)
plot(sil, col=1:5, main="", border=NA)

summary(sil)

## Silhouette of 182 units in 5 clusters from silhouette.default(x = groups, dist = d) :
##  Cluster sizes and average silhouette widths:
##         37         27         30         41         47 
## 0.04886944 0.13933644 0.19057532 0.20516587 0.39670270 
## Individual silhouette widths:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.08574  0.09069  0.19322  0.21068  0.30593  0.54489

Vemos que el average silhouette width no es demasiado alto; probaremos bajando el número de clusters.

Además, nos damos cuenta en este punto de que las transformaciones que hemos realizado previamente de las variables han “estropeado” los datos, en el sentido de que no queda tan clara la división, basada en los dos primeros componentes principales, como cuando utilizábamos los datos sin transformar. Por lo tanto, recogemos cuerda y utilizaremos los datos sin transformar, pues con estos los resultados son bastante más interpretables. Eso puede suceder porque, pese a que las variables por separado se comporten de forma más parecida a la normal, la distribución multidimensional de todo el dataset no lo haga.

X <- data.frame(scale(X_unscaled))

fit = kmeans(X, centers=4, nstart=100)
groups = fit$cluster
#groups
barplot(table(groups), col="blue")

centers=fit$centers
#centers

# clusplot
fviz_cluster(fit, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X, method="euclidean")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

summary(sil)

## Silhouette of 182 units in 4 clusters from silhouette.default(x = groups, dist = d) :
##  Cluster sizes and average silhouette widths:
##        29        36        84        33 
## 0.1839574 0.3919825 0.2296468 0.1496398 
## Individual silhouette widths:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.07903  0.14221  0.23234  0.23997  0.32905  0.54965

Sube un poco el average silhouette width. Seguimos bajando en el número de centroides para ver si mejora.

fit = kmeans(X, centers=3, nstart=100)
groups = fit$cluster
#groups
barplot(table(groups), col="blue")

centers=fit$centers


# clusplot
fviz_cluster(fit, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X, method="euclidean")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

summary(sil)

## Silhouette of 182 units in 3 clusters from silhouette.default(x = groups, dist = d) :
##  Cluster sizes and average silhouette widths:
##        39        55        88 
## 0.3987943 0.3057256 0.2464965 
## Individual silhouette widths:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.0734  0.1934  0.3212  0.2970  0.4099  0.5617

Vemos que con 3 clusters el average silhouette width sigue mejorando; además, estos clusters, visualmente, por los países que forman los clusters, podríamos decir que representarían el primer, segundo y tercer mundo respectivamente, aunque antes de saltar a las conclusiones debemos hacer un estudio detallado de las características de cada cluster para poder afirmar esto.

centers

##   gdp_per_cap  inflation gni_per_cap     urban_p productivity credit_inf
## 1   1.6211935 -0.5446340   1.6543070  1.10617407    1.4874171  0.5075753
## 2  -0.6356268  0.4246104  -0.6539681 -0.87108893   -0.8400000 -0.8068516
## 3  -0.3212167 -0.0240096  -0.3244287  0.05419435   -0.1341962  0.2793341
##     renew_en renew_electr electricity_2016 clean_fuels_2016 inf_mort_rate
## 1 -0.5579749   -0.1685010        0.5993024        0.8531021    -0.9127375
## 2  1.0390537    0.3440539       -1.2998611       -1.3107411     1.2595349
## 3 -0.4021242   -0.1403571        0.5468133        0.4411339    -0.3827006
##   sanitation_services deaths_household_pollution      water  fert_rate
## 1           0.8557941                 -0.9417420  0.7846983 -0.8116898
## 2          -1.3447276                  1.1469979 -1.2734136  1.2176491
## 3           0.4611823                 -0.2995108  0.4481194 -0.4013045
##     life_exp urban_pop_growth    emp_agr    eco_agr vulner_empl
## 1  1.1401976       -0.5241457 -0.9815974 -0.8516165  -1.0470265
## 2 -1.2200850        0.9950285  1.1614203  1.0139171   1.1546237
## 3  0.2572383       -0.3896010 -0.2908616 -0.2562773  -0.2576167
##       acc_own
## 1  1.30387480
## 2 -0.85861219
## 3 -0.04122098

Como vemos, el primer cluster tiene un gdp per capita bastante inferior ala media, una inflación bastante superior, gni per capita muy bajo, población urbana, productividad e información al crédito aún más bajas. Por otro lado, el uso de energías renovables es muy alto (siempre con respecto a la media), el acceso a electricidad renovable también alto aunque no tanto, el uso de electricidad y el acceso a clean fuels en 2016 ambos muy bajos, una alta tasa de mortalidad infantil (altísima), un acceso a los servicios sanitarios bajísimo con respecto a los demás, muertes por household pollution muy altas, bajísimo acceso a agua limpia, una alta tasa de fertilidad, muy baja esperanza de vida, alto crecimiento de la población urbana, alto empleo y economía agrícola, un muy alto número de empleos vulnerables, y muy baja account ownership. Como podemos ver, este primer cluster correspondería al tercer mundo claramente. Serían aquellos países que peor posicionados están con respecto al spread.

Por otro lado, el segundo cluster sigue teniendo un gdp per capita bajo, una inflación cercana a la media, un gni per capita bajo (aunque no tan bajo como el de los países del cluster 1), una poblacion urbana normal (media), productividad cercana a la media, información al crédito alta, aunque no tan alta como en el cluster 3, bajo uso de energías renovables, bajo uso de electricidad renovable aunque cercano a la media. Por otro lado, estos países tienen alto acceso a la electricidad y a comnustibles limpios en 2016, una mortalidad infantil moderadamente baja, un acceso a servicios sanitarios moderadamente alto, bajas muertes por household pollution, moderadamente alto acceso al agua, una tasa de fertilidad baja, aunque no tan baja como la de los países del cluster 3, un crecimiento de la población urbana moderadamente bajo; tanto en empleo agrícola, como en agricultura agrícola, como en empleos vulnerables tienen una puntuación ligeramente inferior a la media, y un porcentaje de personas con cuenta bancaria muy cercano a la media. Este cluster correspondería a los países del segundo mundo, aquellos que están posicionados más o menos en la media con respecto al spread; ni se están aprovechando tanto de las consecuencias del progreso tecnológico como los del tercer clusters ni están padeciendo tanto estos efectos como los del primer cluster.

Por último, el tercer cluster tiene un gdp per capita, gni per capita, población urbana, productividad muy altas; una inflación bastante inferior a la media. Tiene la mayor información crediticia de todos los clusters, un bajo uso de energías renovables, un uso de electricidad renovable muy similar a la del segundo cluster, un acceso a la electricidad alto y acceso a los combustibles limpios similar al del cluster 2. La mortalidad infantil en este cluster es la más baja de todas, mientras que el acceso a servicios sanitarios es bastante alto, las muertes por household pollution es muy baja, el acceso al agua limpia es el más alto de todos los clusters, la tasa de fertilidad es muy baja, mientras que la esperanza de vida es altísima, muy por encima de los otros dos clusters. Tienen el crecimiento de la población urbana más bajo de todos los clusters, y el peso tanto del trabajo como de la economía agrícola es muy bajo en sus sistemas de producción. El porcentaje de empleo vulnerable es muchisimo más bajo que en el resto de clusters, y llama la atención el porcentaje de personas con cuenta bancaria, en cuyo apartado están claramente posicionados muy positivamente, con mucha diferencia con respecto a los otros dos clusters. Este cluster podría ser identificado, por lo tanto, como aquellos países que mejor posicionados están con respecto al efecto del spread; es decr, el primer mundo.

Estamos cómodos con esta composición de clusters, especialmente después de observar cuál es la media de los países de cada cluster en cada una de las variables introducidas, pues nos cuadra con la idea que teníamos al principio de que probablemente encontráramos tres grupos en este aspecto: aquellos que están quedándose claramente atrás debido al progreso tecnológico, los que están más o menos en la media, es decir que no sufren tanto las consecuencias de este progreso pero tampoco se están distanciando positivamente de los demás debido a este; y por último los países que están sacando una ventaja significativa del resto en la calidad de vida, es decir aquellos países que están muy bien posicionados con respecto al spread.

Vamos a probar con 2 centroides:

fit = kmeans(X, centers=2, nstart=100)
groups = fit$cluster
#groups
barplot(table(groups), col="blue")

centers=fit$centers
#centers

# clusplot
fviz_cluster(fit, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X, method="euclidean")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

summary(sil)

## Silhouette of 182 units in 2 clusters from silhouette.default(x = groups, dist = d) :
##  Cluster sizes and average silhouette widths:
##        59       123 
## 0.3541759 0.4486794 
## Individual silhouette widths:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.01397  0.34020  0.45628  0.41804  0.52190  0.61582

Vemos que con 2 clusters parece que incluso mejora con respecto a 3.

centers

##   gdp_per_cap  inflation gni_per_cap    urban_p productivity credit_inf
## 1  -0.6331667  0.3760846  -0.6501957 -0.8903473   -0.8330055 -0.7459775
## 2   0.3037141 -0.1803983   0.3118825  0.4270771    0.3995717  0.3578266
##     renew_en renew_electr electricity_2016 clean_fuels_2016 inf_mort_rate
## 1  0.9920984    0.3636492       -1.1917260       -1.2550132     1.1963608
## 2 -0.4758846   -0.1744333        0.5716409        0.6019982    -0.5738641
##   sanitation_services deaths_household_pollution     water  fert_rate
## 1          -1.2538326                  1.1193596 -1.183160  1.1334424
## 2           0.6014319                 -0.5369286  0.567532 -0.5436837
##     life_exp urban_pop_growth    emp_agr    eco_agr vulner_empl    acc_own
## 1 -1.1542185        0.9343703  1.1373573  0.9968732   1.1343759 -0.8069055
## 2  0.5536495       -0.4481939 -0.5455616 -0.4781750  -0.5441315  0.3870523

Sin embargo, pese a que matemáticamente parecen encajar mejor dos clusters que tres, al observar la media que tiene cada uno de los clusters en las diferentes variables nos encontraríamos más cómodos teniendo 3 grupos, pues por ejemplo en las variables económicas, en este caso, no se observa tanta diferencia (porque agrupa a los países del segundo mundo en el primero o en el tercero sin pertenecer del todo a ninguno de estos dos, lo cual hace que las diferencias en la media de estas variables no sea tan significativa). Mientras que con 3 clusters sí podíamos identificar una gran diferencia del tercer cluster (primer mundo) con respecto a los otros dos en algunas variables como el porcentaje de adultos con cuenta bancaria, o las variables económicas ya mencionadas antes, además de otras variables que miden la calidad de vida como el acceso al agua o a la sanidad, en este caso las diferencias son más sutiles porque los dos clusters están formados por países más dispares entre sí (Brasil y Luxemburgo, por ejemplo, están en el mismo cluster, y mientras que el primero está aún a años luz (aunque tenga mucho potencial) de posicionarse bien con respecto al spread, el segundo podríamos decir que es el país más desarrollado del mundo, o uno de ellos, junto a Singapur, Islandia…).

En caso de utilizar un k-means cluster, por lo tanto, ¿con cuántos nos quedamos?:

fviz_nbclust(X, kmeans, method = 'wss')

fviz_nbclust(X, kmeans, method = 'silhouette')

Ambos gráficos nos vienen a decir lo que acabábamos de ver nosotros probando y viendo el average silhouette width, es decir, que lo óptimo parecen ser 2 clusters. Esto, de todas formas, es sólo de manera matemática, y como dijimos antes nos sentiríamos más cómodos cogiendo 3 clusters.

Recordando lo que hicimos con los datos en el estudio de las variables y el ajuste de PCA y Factor Analysis, también podríamos eliminar algunas variables de nuestro set de datos que no tienen tanta relación con las demás y podrían estar incluyendo ruido en nuestra muestra.

X_reduced <- data.frame(X)

X_reduced[ , c("inflation", "urban_p", "credit_inf", "renew_electr", "urban_pop_growth")] <- NULL

fviz_nbclust(X_reduced, kmeans, method = 'wss')

fviz_nbclust(X_reduced, kmeans, method = 'silhouette')

Con estos datos reducidos, cogiendo 3 centroides, quedaría asi:

fit = kmeans(X_reduced, centers=3, nstart=100)
groups = fit$cluster
#groups
barplot(table(groups), col="blue")

centers=fit$centers
#centers

# clusplot
fviz_cluster(fit, data = X_reduced, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X_reduced, method="euclidean")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

summary(sil)

## Silhouette of 182 units in 3 clusters from silhouette.default(x = groups, dist = d) :
##  Cluster sizes and average silhouette widths:
##        53        91        38 
## 0.3851991 0.3602663 0.4495979 
## Individual silhouette widths:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.0574  0.2872  0.4304  0.3862  0.5154  0.6165

i=1  # plottinng the centers in cluster 1
barplot(centers[i,], las=2, col="darkblue", main=paste("Cluster", i))

i=2  # plottinng the centers in cluster 1
barplot(centers[i,], las=2, col="darkblue", main=paste("Cluster", i))

i=3

barplot(centers[i,], las=2, col="darkblue", main=paste("Cluster", i))

Sin embargo, pese a observar que los países que forman cada cluster tendrían más en común entre sí, pues el average silhouette width sube, decidimos quedarnos con los datos originales, es decir con la matriz X. Es evidente que al reducir la dimensionalidad de esta matriz (quitando columnas), tendríamos mejores resultados, pero esto puede deberse únicamente al bajo número de observaciones que tenemos, y no tanto a la poca relevancia de esas variables. Por este motivo decidimos trabajar con la matriz completa aunque esto implique incluir un poco de ruido en nuestros datos, lo cual, como comentamos, es inevitable debido a que \(\frac{p}{n}\) es grande.

Using Mahalanobis Distance

Hasta este punto hemos estado haciendo cluster con las distancias euclídeas, pero vamos a probar a hacerlo con la distancia de Mahalanobis.

# k-means with Mahalanobis distance

S_x <- cov(X_unscaled)
iS <- solve(S_x) #para coger la inversa de la matrix
e <- eigen(iS)
V <- e$vectors
B <- V %*% diag(sqrt(e$values)) %*% t(V)
Xtil <- scale(X_unscaled,scale = FALSE)
XS <- Xtil %*% B

# kmeans with 2 clusters
fit = kmeans(XS, centers=2, nstart=100)
groups = fit$cluster
barplot(table(groups), col="blue")

centers=fit$centers
colnames(centers)=colnames(X)
#centers

# clusplot
fviz_cluster(fit, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

fviz_nbclust(XS, kmeans, method = 'wss')

fviz_nbclust(XS, kmeans, method = 'silhouette')

En los gráficos de arriba queda suficientemente claro que realizar un k-means clusters utilizando la distancia de Mahalanobis nos da peores resultados que utilizando la distancia euclídea. Por este motivo no entraremos a explicar las características de cada uno de estos clusters, e iremos directamente a implementar el siguiente método de aprendizaje no supervisado: PAM; no sin antes visualizar los resultados del modelo de cluster por el método de k-means que más nos cuadra de los probados arriba.

if(!require(rworldmap)) {
  
  install.packages("rworldmap")
  
  
}

fit = kmeans(X, centers=3, nstart=100)
groups = fit$cluster

matched <- joinCountryData2Map(data.frame(country=countries, value=groups), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World by K-Means", 
               catMethod = "categorical", colourPalette = "topo")

Como vemos en el mapa del mundo, nos cuadran bastante los países que pertenecen al primer mundo (cluster 3), segundo mundo (cluster 2) y tercer mundo (cluster 1). Probaremos ahora con el siguiente método: PAM.

4. PAM

Using Euclidean Distance

# number of clusters?
fviz_nbclust(X, pam, method = 'wss')

fviz_nbclust(X, pam, method = 'silhouette')

# PAM clustering, based on medioids
fit.pam <- pam(X, k=2, metric = "euclidean") # 2 cluster solution
fit.pam$medoids

##    gdp_per_cap  inflation gni_per_cap    urban_p productivity credit_inf
## 33  -0.6435012 -0.7046321  -0.6525326 -0.1063936   -0.8740631 -1.2085410
## 42   0.3844251 -0.3483195   0.3373432  0.6971845    0.6696516  0.7408338
##      renew_en renew_electr electricity_2016 clean_fuels_2016 inf_mort_rate
## 33  1.5459515    1.3198344       -0.7815197       -1.1109502     1.6683661
## 42 -0.6093636   -0.6762417        0.6553683        0.8557089    -0.9702128
##    sanitation_services deaths_household_pollution      water  fert_rate
## 33          -1.1885475                  1.3631735 -1.1869961  1.4674467
## 42           0.8566682                 -0.8092583  0.7934034 -0.9611414
##      life_exp urban_pop_growth   emp_agr    eco_agr vulner_empl    acc_own
## 33 -1.8001742         0.895438  1.535388  0.8394029   1.3416915 -0.9453004
## 42  0.8915702        -1.048183 -1.014676 -0.7953806  -0.8615597  0.8047001

#fit.pam$clustering


fviz_cluster(fit.pam, data = X, geom = c("point"), ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

Vemos que de nuevo, matemáticamente, parece que lo que más sentido tiene es tomar 2 clusters. Probaremos con otras distancias, así como con 3 clusters, y explicaremos los resultados de aquella composición que más nos encaje para el objetivo de este estudio.

# PAM clustering, based on medioids
fit.pam <- pam(X, k=3, metric = "euclidean") # 2 cluster solution
fit.pam$medoids

##    gdp_per_cap  inflation gni_per_cap    urban_p productivity credit_inf
## 33  -0.6435012 -0.7046321  -0.6525326 -0.1063936   -0.8740631 -1.2085410
## 38  -0.5477250 -0.6763866  -0.5564911  0.3731611   -0.6336868  0.4159380
## 51   0.8075529 -0.4456493   0.8688122  0.9823251    1.1591846  0.7408338
##      renew_en renew_electr electricity_2016 clean_fuels_2016 inf_mort_rate
## 33  1.5459515   1.31983443       -0.7815197       -1.1109502     1.6683661
## 38 -0.1989357  -0.40468249        0.3894632        0.1634364    -0.3470056
## 51 -0.5595744   0.05001927        0.6553683        0.9318907    -0.9702128
##    sanitation_services deaths_household_pollution        water  fert_rate
## 33          -1.1885475                 1.36317353 -1.186996132  1.4674467
## 38          -0.3022874                -0.01661423  0.001243579 -0.2672591
## 51           0.8907551                -0.98540139  0.793403387 -1.0999179
##      life_exp urban_pop_growth     emp_agr    eco_agr vulner_empl
## 33 -1.8001742       0.89543799  1.53538810  0.8394029  1.34169151
## 38  0.1366442      -0.04692393  0.04749208 -0.2288715  0.02836333
## 51  1.4237312      -0.81259299 -0.95447259 -0.7872876 -0.94770847
##       acc_own
## 33 -0.9453004
## 38 -0.5864026
## 51  1.2860648

#fit.pam$clustering

groups <- fit.pam$cluster

fviz_cluster(fit.pam, data = X, geom = c("point"), ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X, method="euclidean")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

Using Manhattan Distance

dis <- dist(X, method = "manhattan")
# number of clusters?
fviz_nbclust(X, pam, method = 'wss', diss=dis)

fviz_nbclust(X, pam, method = 'silhouette', diss=dis)

fit.pam <- pam(X, k=2, metric = "manhattan") # 2 cluster solution
#fit.pam$medoids
#fit.pam$clustering

groups <- fit.pam$cluster

fviz_cluster(fit.pam, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X, method="manhattan")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

Vemos que en el caso de PAM utilizando la distancia de Manhattan, lo óptimo (matemáticamente) también sería coger dos clusters. Con la distancia de Manhattan el average silhouette width sale algo más alto que con la distancia euclídea, sin importar el número de clusters (para todos en general es más alta).

fit.pam <- pam(X, k=3, metric = "manhattan") 

groups <- fit.pam$cluster

fviz_cluster(fit.pam, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X, method="manhattan")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

Con 3 clusters nos ocurre lo mismo que comentábamos antes, y es que nos encajan más los resultados que cuando agrupamos a lo países solamente en 2 grupos. Vemos que utilizando el método de PAM con la distancia Manhattan mejoramos los resultados para este cluster de 3; midiéndolo como el average silhouette width y comparándolo con PAM cuando usábamos las distancias de Manhattan y con k-means usando las distancias euclídea y de Mahalanobis.

Using Mahalanobis Distance

S_x <- cov(X_unscaled)
iS <- solve(S_x) #para coger la inversa de la matrix
e <- eigen(iS)
V <- e$vectors
B <- V %*% diag(sqrt(e$values)) %*% t(V)
Xtil <- scale(X_unscaled,scale = FALSE)
XS <- Xtil %*% B


fit.pam <- pam(XS, k = 3, metric = "euclidean")

groups <- fit.pam$cluster

fviz_cluster(fit.pam, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

d <- dist(X, method="euclidean")  
sil = silhouette(groups, d)
plot(sil, col=1:4, main="", border=NA)

De nuevo utilizar la distancia de Mahalanobis parece una idea pésima, pues no se distinguen bien los grupos, no agrupa tan bien como el método de PAM usando las distancias Euclidea y de Manhattan. El mejor resultado de cluster para la técnica PAM, es por tanto cogiendo 3 clusters (matemáticamente serían 2, pero por nuestro objetivo del trabajo, como dijimos antes, nos encajan mejor 3), con la distancia de Manhattan.

Vamos a volver a visualizar estos resultados y a explicar cuáles son las caractrísticas de cada uno de estos clusters.

fit.pam <- pam(X, k=3, metric = "manhattan") 
fit.pam$medoids

##    gdp_per_cap   inflation gni_per_cap     urban_p productivity credit_inf
## 33  -0.6435012 -0.70463209  -0.6525326 -0.10639363   -0.8740631 -1.2085410
## 83  -0.4445131  0.03126667  -0.4527898 -0.08911238   -0.5871169  0.7408338
## 56   1.3681889 -0.62749405   1.5045717  0.97800480    1.4412508  0.4159380
##      renew_en renew_electr electricity_2016 clean_fuels_2016 inf_mort_rate
## 33  1.5459515    1.3198344       -0.7815197       -1.1109502      1.668366
## 83 -0.5414891   -0.7114842        0.5907409        0.6797209     -0.442497
## 56 -0.6557915   -0.5388849        0.6553683        0.9318907     -0.924980
##    sanitation_services deaths_household_pollution      water  fert_rate
## 33          -1.1885475                  1.3631735 -1.1869961  1.4674467
## 83           0.3794512                 -0.6918295  0.3973235 -0.5448120
## 56           0.8566682                 -0.9854014  0.7934034 -0.6142002
##      life_exp urban_pop_growth    emp_agr    eco_agr vulner_empl
## 33 -1.8001742        0.8954380  1.5353881  0.8394029  1.34169151
## 83  0.5636106       -0.6947977 -0.3524337 -0.4231032  0.03917959
## 56  1.3556641       -0.7536954 -1.0060759 -0.8601245 -1.12229910
##       acc_own
## 33 -0.9453004
## 83 -0.8579625
## 56  1.2954307

fviz_cluster(fit.pam, data = X, geom = c("point"),ellipse.type = 'norm', pointsize=1)+
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

Vemos que el primer cluster tiene un gdp per capita muy bajo (el más bajo de los 3 clusters), una inflación muy baja también, gni per capita muy bajo, población urbana baja pero cercana a la media, productividad muy baja (la más baja de todas), información al crédito muy muy baja, un altísimo uso de energías y de electricidad renovable, sin embargo bajísimo acceso a la electricidad y a combustibles limpios en el año 2016, una tasa de mortalidad infantil altísima, servicios sanitarios bajísimos, muy inferiores a la media, muertes por household pollution muy altas, bajíismo acceso al agua, altísima tasa de fertilidad, bajísima esperanza de vida en comparación con lo otros dos clusters; además, es el cluster que tiene una media de crecimiento de la población urbana más alta, con un altísimo peso del empleo y economía agrícola, un gran número de empleos vulnerables y un bajísimo porcentaje medio de account ownership. Vemos por lo tanto que este primer cluster representaría el tercer mundo, o aquellos países que peor posicionados están con respecto al efecto del spread. El cluster 2 vemos que tiene también unas variables económicas bajas, aunque no tan bajas como las del cluster 1… Muy similar al análisis que realizamos antes para el 3-mean custering con distancia euclídea; por lo tanto podríamos decir que este segundo cluster está formado por los países del segundo mundo, o aquellos países que no son por tanto ni los mejor posicionados con respecto al spread ni tampoco se están quedando atrás debido a este. EL último y tercer cluster está formado por los países que mejor puntuación media tienen en las variables económicas, así como en calidad de vida y sanidad (agua, servicios sanitarios…), como ya hemos dicho antes estos resultados son casi iguales que los que teníamos con el método de 3-mean clustering usando distancias euclídeas. Podemos identificar este tercer cluster con los países del primer mundo, o aquellos que estan mejor posicionados con respecto al spread.

Antes de pasar al siguiente método de clustering, vamos a dibujar en el mapa del mundo a los países según el grupo al que pertenecen por el mejor modelo de PAM clustering; es decir, 3 clusters utilizando la distancia de Manhattan.

fit.pam <- pam(X, k=3, metric = "manhattan") 

groups <- fit.pam$cluster

matched <- joinCountryData2Map(data.frame(country=countries, value=groups), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using PAM", 
               catMethod = "categorical", colourPalette = "topo")

Vemos que por el método PAM los clusters salen muy parecidos, excepto alguna salvedad como Chile, que pasa de estar en el segundo mundo al primero.

5. Model-Based Clustering (Gaussian Mixtures)

Primero vamos a dejar que sea el propio modelo el que optimice el número de clusters (anticipamos que seguramente lo divida en 2, como hemos visto antes parece que al dividirlos en 2 matemáticamente es con el número de clusters que más tienen en común los países de cada cluster, teóricamente); luego después ajustaremos un modelo de 3 clusters, por los motivos que ya explicamos arriba, pues es con el número de grupos que más cómodos nos sentimos.

require(mclust)
fit.m <- Mclust(X)
summary(fit.m) # seems 2 clusters fit well

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVE (ellipsoidal, equal orientation) model with 2 components: 
## 
##  log.likelihood   n  df       BIC       ICL
##       -1952.266 182 295 -5439.714 -5440.387
## 
## Clustering table:
##   1   2 
## 106  76

# we can also specify tghe number of groups
# fit.m <- Mclust(X, G=3)

plot(fit.m, what = "classification", main="")

plot(fit.m, what = "BIC", main="")

# By default, the models considered are:
# "EII": spherical, equal volume 
# "VII": spherical, unequal volume 
# "EEI": diagonal, equal volume and shape
# "VEI": diagonal, varying volume, equal shape
# "EVI": diagonal, equal volume, varying shape 
# "VVI": diagonal, varying volume and shape 
# "EEE": ellipsoidal, equal volume, shape, and orientation 
# "EEV": ellipsoidal, equal volume and equal shape
# "VEV": ellipsoidal, equal shape 
# "VVV": ellipsoidal, varying volume, shape, and orientation

groups = fit.m$classification

centers = fit.m$parameters$mean

i=1  # plottinng the centers in cluster 1
barplot(centers[,i], las=2, col="darkblue", main=paste("Cluster", i))

# insights?

i=2  # plottinng the centers in cluster 2
barplot(centers[,i], las=2, col="darkblue", main=paste("Cluster", i))

# scatter plot:
fviz_cluster(fit.m,geom = c("point"), ellipse.type = 'norm',pointsize=1) +
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

Ahora vamos a ajustar un modelo de 3 clusters:

fit.m <- Mclust(X, G=3)
summary(fit.m)

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEV (ellipsoidal, equal shape) model with 3 components: 
## 
##  log.likelihood   n  df      BIC       ICL
##       -1184.192 182 718 -6104.86 -6104.961
## 
## Clustering table:
##  1  2  3 
## 73 39 70

plot(fit.m, what = "classification", main="")

plot(fit.m, what = "BIC", main="")

groups = fit.m$classification

centers = fit.m$parameters$mean

i=1  # plottinng the centers in cluster 1
barplot(centers[,i], las=2, col="darkblue", main=paste("Cluster", i))

# insights?

i=2  # plottinng the centers in cluster 2
barplot(centers[,i], las=2, col="darkblue", main=paste("Cluster", i))

# insights?

i = 3
barplot(centers[,i], las=2, col="darkblue", main=paste("Cluster", i))

# scatter plot:
fviz_cluster(fit.m,geom = c("point"), ellipse.type = 'norm',pointsize=1) +
  theme_minimal()+geom_text(label=countries,hjust=0, vjust=0,size=2,check_overlap = T)+scale_fill_brewer(palette="Paired")

Analizaremos ahora las características de los clusters basándonos en los 3 gráficos de barras que hemos dibujado arriba. En ellos se ve la media (recordemos que todas las variables están normalizadas) que tiene cada uno de los clusters en las 21 variables sobre los países que estamos considerando.

El primer cluster vemos que tiene unas variables económicas (gdp per capita, gni per capital, productividad, información crediticia…) bastante bajas; por debajo de la media. Llaman la atención el bajísimo acceso a la electricidad y a combustibles limpios en el año 2016, la alta tasa de mortalidad infantil, el bajo acceso a servicios sanitarios y al agua, la alta tasa de fertilidad… En general, vemos de nuevo que este primer cluster correspondería a aquellos países que pertenecen al tercer mundo, es decir lo forman los países peor posicionados con respecto al spread.

El segundo cluster está claro que pertenece a los países que forman parte del primer mundo, es decir aquellos países que están mejor posicionados con respecto al spread. Esto lo vemos claro en el gráfico de barras, pues es el grupo que mayor media tiene en las variables económicas y otras variables que miden la calidad de vida, mientras que tiene una tasa bajísima de mortalidad infantil, muy pocos empleos vulnerables etc.

Por último, el tercer y último cluster está entre medias de los otros dos en casi todas las variables, podemos ver que no puntúa tan alto en las variables “buenas” (aquellas asociadas con buen posicionamiento con respecto al spread) como el segundo cluster, ni tan mal como el primero; por otro lado, no puntúa tan bajo en las variables “malas” como el segundo, ni tan alto como el primero. Por este motivo esto nos lleva a decir que este tercer cluster son los países que están posicionados de forma intermedia con respecto al spread.

Ahora visualizamos los resultados de este último modelo de cluster en el mapa del mundo, para ver si existen diferencias significativas con respecto al mejor modelo de PAM y de k-means. Primero cogeremos los grupos, y cambiaremos el orden de estos para que visualmente sea más fácilmente comparable con los otros dos gráficos que dibujamos antes, pues recordemos que en este caso el primer mundo por ejemplo en lugar de ser el cluster 3 es el 2.

matched <- joinCountryData2Map(data.frame(country=countries, value=new_groups), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using PAM", 
               catMethod = "categorical", colourPalette = "topo")

Vemos que en este caso los resultados también se asemejan a los obtenidos por los otros dos métodos de clustering; se asemeja más a la separación de grupos que había hecho el k-means. Sin embargo, este modelo es algo más fino y países como India los mete en el cluster 1, es decir en los países del tercer mundo.

Por último, antes de comparar todos los métodos de cluster y de quedarnos con el modelo que más nos convenza, vamos a implementar Hierarchical Clustering.

6. Hierarchical Clustering.

En el caso del hierarchical clustering, y dado que ya hemos comentado sobradamente que 3 clusters nos encajan mejor que 2 para el objetivo de este estudio, vamos a probar solamente haciendo 3 clusters (cortando el árbol de tal forma que dejemos 3 grupos, ya que en este método de cluster, en función de donde “cortes el árbol” obtienes un número u otro de clusters).

Hierarchical Clustering with Euclidean Distance

if(!require(ape)) {
  
  install.packages("ape")
  
  require(ape)
  
}

## Loading required package: ape

## Warning: package 'ape' was built under R version 3.5.2

d <- dist(X, method = "euclidean") # distance matrix using euclidean distance
# we can choose any of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski" 


fit.h <- hclust(d, method="ward.D2") # cluster using Ward linkage (based on minimizing variance)

dendro <- as.dendrogram(fit.h)

fviz_dend(fit.h, k=3, rect = TRUE) # dendrogam

groups.h <- cutree(fit.h, k=3) 

X_grouped <- X

X_grouped$group <- groups.h

En este caso dibujamos el mapa del mundo con los clusters para cada modelo de Hierarchical Clustering, pues es la forma más clara de determinar, visualmente, si nos encajan los resultados obtenidos en cada caso.

matched <- joinCountryData2Map(data.frame(country=countries, value=groups.h), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

Vemos que los resultados son muy similares a los obtenidos por el método de k-means y PAM anteriormente (más parecido al primero).

Haremos la media de las variables para cada uno de los clusters para poder determinar con algo más de claridad las características de cada uno de ellos.

require(dplyr)

X_grouped %>%
  
  group_by(group) %>%
  
  summarise_all(funs(mean))

## # A tibble: 3 x 22
##   group gdp_per_cap inflation gni_per_cap urban_p productivity credit_inf
##   <int>       <dbl>     <dbl>       <dbl>   <dbl>        <dbl>      <dbl>
## 1     1      -0.634    0.405       -0.652  -0.873      -0.832      -0.758
## 2     2      -0.253   -0.0575      -0.255   0.110      -0.0769      0.266
## 3     3       1.86    -0.555        1.90    1.24        1.71        0.578
## # ... with 15 more variables: renew_en <dbl>, renew_electr <dbl>,
## #   electricity_2016 <dbl>, clean_fuels_2016 <dbl>, inf_mort_rate <dbl>,
## #   sanitation_services <dbl>, deaths_household_pollution <dbl>,
## #   water <dbl>, fert_rate <dbl>, life_exp <dbl>, urban_pop_growth <dbl>,
## #   emp_agr <dbl>, eco_agr <dbl>, vulner_empl <dbl>, acc_own <dbl>

Vemos que el primer grupo son los países del tercer mundo, con un gdp per capita y un gni per capita muy bajos, con poquísima población urbana (siempre hablamos en medias y con respecto a los otros grupos), con bastante inflación (mucha más que los otros dos grupos); además es el grupo con menor información de crédito, y menor acceso a electricidad y a combustibles limpios en 2016. Curiosamente los países del tercer mundo son los que mayor uso de electricidad y energías renovables tienen. Sin embargo, tienen una mortalidad infantil altísima, unos servicios sanitarios muy bajos, así como bajísimo acceso al agua, muchas muertes por household pollution, una esperanza de vida muy inferior a la de los otros grupos, un gran crecimiento de la población urbana… El empleo y la economía agrícola tienen un gran peso en las sociedades que forman parte de este primer cluster, que tienen además un gran número de empleos vulnerables y un bajo porcentaje de adultos con cuenta bancaria en su posesión.

El segundo grupo serían los países que forman el segundo mundo, ya que vemos que tienen las características “buenas” (sanidad, agua, variables económicas, posesión de cuenta bancaria, esperanza de vida) más altas que el cluster 1 pero más bajas que el cluster 2. Además, las características “negativas” (muertes por contaminación, empleo y economía agrícolas, empleo vulnerable, inflación…) son más bajas en este segundo cluster que en el primero, pero más altas que en el tercero.

Por último, el tercer y último cluster lo forman aquellos países que tienen mayor puntuación en las variables de índole positiva (productividad, esperanza de vida… todas las mencionadas antes), mientras que tienen la puntuación más baja en las variables de índole negativa, como mortalidad infantil o empleo vulnerable (además de las ya mencionadas previamente como negativas). Podemos concluir por tanto que este tercer cluster lo forman los países del primer mundo, es decir aquellos que mejor posicionados están con respecto al spread.

Ahora probaremos cogiendo otras distancias, en este caso la distancia de Manhattan.

Hierarchical Clustering Using Manhattan Distance

d <- dist(X, method = "manhattan") # distance matrix using euclidean distance
# we can choose any of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski" 

fit.h <- hclust(d, method="ward.D2")

dendro <- as.dendrogram(fit.h)



fviz_dend(fit.h, k=3, rect = TRUE)

groups.h <- cutree(fit.h, k=3)

matched <- joinCountryData2Map(data.frame(country=countries, value=groups.h), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

Utilizando las distancias de Manhattan nos convence menos el hierarchical cluster, pues como vemos algunos países como Rusia o Argentina (el uno con un régimen político de dudosa legitimidad moral y el otro con problemas económicos cíclicos) caerían dentro de los países del primer mundo, lo cual no nos cuadra. Probaremos ahora con otro tipo de distancia.

Hierarchical Clustering Using Maximum Distance.

d <- dist(X, method = "maximum") # distance matrix using euclidean distance
# we can choose any of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski" 

# With the input, we select the linkage for the cluster joining (simple, complete, average, Ward)
fit.h <- hclust(d, method="ward.D2") # cluster using Ward linkage (based on minimizing variance)

dendro <- as.dendrogram(fit.h)

fviz_dend(fit.h, k=3, rect = TRUE) # dendrogam

groups.h <- cutree(fit.h, k=3) 

matched <- joinCountryData2Map(data.frame(country=countries, value=groups.h), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

Vemos que en este caso, utilizando como medida de distancias “Maximum”, separa los países parecido a los otros métodos utilizados, con algunas salvedades como que más países africanos caerían dentro del segundo mundo, y más países del primer mundo también quedarían dentro de este segundo mundo.

Hierarchical Clustering Using Minkowski Distance

d <- dist(X, method = "minkowski") # distance matrix using euclidean distance
# we can choose any of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski" 

# With the input, we select the linkage for the cluster joining (simple, complete, average, Ward)
fit.h <- hclust(d, method="ward.D2") # cluster using Ward linkage (based on minimizing variance)

dendro <- as.dendrogram(fit.h)

# Then, we can visualize the dendogram
#plot(fit.h, labels = countries, hang = -1) 

# And we can also cut the tree in our desired number of groups
#rect.hclust(fit.h, k=3, border="red")

# Nicer but slower
fviz_dend(fit.h, k=3, rect = TRUE) # dendrogam

# Instead of visualizing, we can create the groups by choosing our convenient number of groups 
groups.h <- cutree(fit.h, k=3) 

#fviz_silhouette(fit.h)

matched <- joinCountryData2Map(data.frame(country=countries, value=groups.h), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

Vemos que apenas cambia nada (salvo algunos países como India) con respecto al Hierarchical Clustering utilizando las distancias máximas.

7(bis). Advanced Hierarchical Clustering.

Using Euclidean Distance.

require(factoextra)

fit.hc <- eclust(X, 'hclust', k=3, hc_metric = 'euclidean', hc_method = "ward.D2", graph=F)
groups = fit.hc$cluster
fviz_dend(fit.hc, rect = TRUE) # dendrogam

fviz_silhouette(fit.hc)

##   cluster size ave.sil.width
## 1       1   57          0.30
## 2       2   93          0.27
## 3       3   32          0.44

matched <- joinCountryData2Map(data.frame(country=countries, value=groups), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

Los resultados de este cluster nos salen muy parecidos a los que obtuvimos antes, como podemos ver en el mapa del mundo.

Using Manhattan Distance.

fit.hc <- eclust(X, 'hclust', k=3, hc_metric = 'manhattan', hc_method = "ward.D2", graph=F)
groups = fit.hc$cluster
fviz_dend(fit.hc, rect = TRUE) # dendrogam

fviz_silhouette(fit.hc)

##   cluster size ave.sil.width
## 1       1   56          0.34
## 2       2   65          0.26
## 3       3   61          0.37

matched <- joinCountryData2Map(data.frame(country=countries, value=groups), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

El average silhouette width sale algo mejor que en el anterior caso, aunque solo 0.01 más; sin embargo vemos en el mapa que nos cuadra menos con la idea que tenemos, pues pasa algo parecido a lo que vimos antes, y es que países como Rusia y Argentina entran dentro del cluster 3 (primer mundo), lo cual no es muy consistente con las ideas que tenemos de estos países.

Using Canberra Distance.

fit.hc <- eclust(X, 'hclust', k=3, hc_metric = 'canberra', hc_method = "ward.D2", graph=F)
groups = fit.hc$cluster
fviz_dend(fit.hc, rect = TRUE) # dendrogam

fviz_silhouette(fit.hc)

##   cluster size ave.sil.width
## 1       1   49          0.47
## 2       2   49         -0.03
## 3       3   84          0.41

matched <- joinCountryData2Map(data.frame(country=countries, value=groups), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

Utilizando la distancia de Canberra los resultados son aún menos coherentes, pues muchos países que según este modelo de clustering entrarían en el primer mundo (como Brasil), está claro que no pueden formar parte del mismo grupo que países como Islandia, Singapur, Luxemburgo…

Using the “Complete” Method instead of Ward.D2.

Ahora vamos a probar, ya como último modelo de clustering (antes de comparar matemáticamente, con algunas estadísticas, todos los modelos de cluster que nos han convencido), el hierarchical clustering utilizando distancias euclídeas pero con el método de Complete.

fit.hc <- eclust(X, 'hclust', k=3, hc_metric = 'euclidean', hc_method = "complete", graph=F)
groups = fit.hc$cluster
fviz_dend(fit.hc, rect = TRUE) # dendrogam

fviz_silhouette(fit.hc)

##   cluster size ave.sil.width
## 1       1   65          0.28
## 2       2  109          0.28
## 3       3    8          0.43

matched <- joinCountryData2Map(data.frame(country=countries, value=groups), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched, nameColumnToPlot="value", mapTitle="Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

Vemos que con este método se quedan solos en el cluster 3 (primer mundo) los países con variables económicas muy por encima del resto, teniendo menos peso a la hora de formar parte de ese cluster el resto de variables del nivel de vida que comentamos antes.

7. Comparing The Best Clustering Models.

Como para sacar estas estadísticas necesitamos utilizar una única matriz de distancias para todos los clusters, y la matriz de distancias que mejor resultados nos ha dado en general ha sido utilizando la distancia euclídea (excepto en el caso de PAM que nos cuadraban mejor los resultados con la distancia de Manhattan), será esta la distancia utilizada para sacar la matriz de distancias.

Con la función cluster.stats comparamos el acuerdo que hay entre 2 composiciones de cluster. Con el corrected Rand Index medimos este nivel de acuerdo, de -1 (no acuerdo en absoluto) a 1 (totalmente de acuerdo). Por otro lado, cuanto más cerca esté el Meila´s VI del 0, más similares son los clusters.

fit.k_means = kmeans(X, centers=3, nstart=100)

groups_kmeans <- fit.k_means$cluster

fit.pam <- pam(X, k=3, metric = "euclidean")

groups_pam <- fit.pam$clustering

fit.m_clust <- Mclust(X, G=3)

groups_mclust <- fit.m$classification

d <- dist(X, method = "euclidean")

fit.h <- hclust(d, method="ward.D2")

groups.h <- cutree(fit.h, k=3)

Comparación K-Means vs PAM

if(!require(fpc)) {
  
  install.packages("fpc")
  
  library(fpc)
}

## Loading required package: fpc

clust_stats = cluster.stats(d, groups_kmeans, groups_pam)

# Among the values returned by previous function, there are two indexes to assess the similarity of two clustering:
#      Corrected Rand index and VI Score
# Rand Index should be maximized and VI score should be minimized
clust_stats$corrected.rand

## [1] 0.5737718

clust_stats$vi

## [1] 0.7556367

Para poder poner en perspectiva los resultados mostrados arriba, es importante comparar todo el resto de clusters con todos.

Comparación K-Means vs Hierarchical Clustering

clust_stats <- cluster.stats(d, groups_kmeans, groups.h)

clust_stats$corrected.rand

## [1] 0.8568725

clust_stats$vi

## [1] 0.3373642

Vemos que el K-Means y el Hierarchical Clustering se parecen mucho, pues tienen un Corrected Rand Index de 0.86, bastante por encima del 0.57 del K-means y PAM; el Meila´s VI es bastante inferior a este anterior caso. Podemos decir, por tanto, que el K-Means y el Hierarchical Clustering tienen bastante similaridad, agrupan a la mayoría de los países en el mismo cluster.

Comparación PAM vs Hierarchical Clustering

clust_stats <- cluster.stats(d, groups_pam, groups.h)

clust_stats$corrected.rand

## [1] 0.4903959

clust_stats$vi

## [1] 0.8554803

Sin embargo, PAM y Hierarchical Clustering no se parecen demasiado, menos que K-Means y PAM y mucho menos que K-Means y Hierarchical Clustering. De momento parece que entre K-Means y PAM nos quedaríamos con el primero; los resultados de ambos debieran ser muy parecidos, con la diferencia de que PAM encuentra sus medioides entre los países que existen, mientras que los centroides de K-Means son ficticios (no son observaciones reales); esto puede significar una mejor interpretabilidad de los clusters de PAM, sin embargo puede llevar (especialmente al tratar con datos de países, pues hay pocas observaciones y los medioides no tienen tanta libertad para elegirse) a una peor agrupación de los países.

Comparación Model-Based Clustering vs. Hierarchical Clustering.

clust_stats <- cluster.stats(d, groups_mclust, groups.h)

clust_stats$corrected.rand

## [1] 0.5890726

clust_stats$vi

## [1] 0.8039513

Vemos que estos dos tipos de cluster se parecen algo más que PAM y Hierarchical Clustering y que PAM y K-Means, sin embargo se parecen bastante menos que Hierarchical Clustering y K-Means.

COmparación K-Means vs. Model-Based Clustering.

clust_stats <- cluster.stats(d, groups_kmeans, groups_mclust)

clust_stats$corrected.rand

## [1] 0.5576133

clust_stats$vi

## [1] 0.913482

Vemos que la similaridad de estos dos clusters es muy similar a la de K-Means y PAM; es bastante inferior a la de K-Means y Hierarchical Clustering.

Comparación PAM vs. Model-Based Clustering

clust_stats <- cluster.stats(d, groups_pam, groups_mclust)

clust_stats$corrected.rand

## [1] 0.4140492

clust_stats$vi

## [1] 1.042174

Como vemos, PAM no tiene mucho que ver con Model-Based Clustering, en términos de los resultados de la agrupación, pues estas dos agrupaciones son bastante diferentes, a juzgar por el corrected rand index y el Meila´s VI score. Podemos determinar que PAM no tiene demasiado en común con ninguno de los otros métodos de Clustering.

8. Conclusiones

Finalmente, se concluye que los dos modelos de Clustering que más tienen en común son K-Means y Hierarchical Clustering; serán por lo tanto los dos modelos de clustering con los que nos quedamos. Volveremos a graficar sus respectivos resultados sobre el mapa del mundo para visualizar las diferencias:

matched1 <- joinCountryData2Map(data.frame(country=countries, value=groups.h), joinCode="NAME", nameJoinColumn="country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched1, nameColumnToPlot="value", mapTitle="Final Clusters for Countries in the World using Hierarchical Clustering", 
               catMethod = "categorical", colourPalette = "topo")

matched2 <- joinCountryData2Map(data.frame(country = countries, value = groups_kmeans), joinCode = "NAME",
                                nameJoinColumn = "country")

## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your data

mapCountryData(matched2, nameColumnToPlot="value", mapTitle="Final Clusters for Countries in the World using K-Means", 
               catMethod = "categorical", colourPalette = "topo")

Vemos que las diferencias entre ambos métodos de cluster son poco significativas. Algunas que nos llaman la atención por ejemplo son:

Portugal: está en el segundo mundo (cluster 2) según el Hierarchical Clustering y en el primer mundo (cluster 3) según el K-Means.
Grecia: caso similar al de Portugal.
India: está en el tercer mundo (cluster 1) según el Hierarchical Clustering y en el segundo mundo (cluster 2) según el K-Means.
Eslovenia, Estonia y Grecia: están en el primer mundo (cluster 3) según el K-Means y en el segundo mundo (cluster 2) según el Hierarchical Clustering.

Vemos, por tanto, que pese a ser similares, el Hierarchical Clustering es algo más exigente a la hora de incluir un país en un cluster avanzado, y tiende a meter a los países en clusters más bajos que el K-Means. En las discrepancias que vemos, nos decidimos quedar con el Hierarchical Clustering, pues nos sentimos más cómodos incluyendo a países como Grecia o Eslovenia en el cluster 2 que en el cluster 1 (no son los países que mejor posicionados están con respecto al spread, y tienen gran diferencia, especialmente en los últimos años – piense el lector en el caso de Grecia) con respecto a los países más desarrollados como Islandia, Singapur o Luxemburgo. Esto es coherente pues generalmente el Hierarchical Clustering es algo más fino que el K-Means clustering; es computacionalmente más costoso, pero esto con los datos que tratamos en este caso no es un problema, pues el número de países existentes es relativamente limitado. En algunos casos nos surgen dudas, como es el caso de Estonia, pues es un país que en los últimos años ha llevado a cabo reformas de transformación digital; además es el país con mayor número de puestos de carga de coches eléctricos de Europa. Sin embargo, históricamente no ha tenido una calidad de vida como la que han gozado los habitantes de países como Luxemburgo. En cualquier caso, los indicadores de calidad de vida como fertility rate o esperanza de vida, así como el acceso a una sanidad de calidad, tardan varios años en estabilizarse, y transformaciones radicales como la observada en Singapur son sucesos remotos. Por este motivo, consideramos que de momento Estonia pertenece a los países del segundo mundo, aquellos que están bien posicionados con respecto al spread en comparación con los países menos desarrollados (cluster 1 – ver gráficos arriba), pero a los que aún les falta para llegar a posicionarse como los países del primer mundo (cluster 3).

Las descripciones de las características de cada cluster (de los dos de arriba) están descritas en sus correspondientes apartados.

IMPLEMENTING UNSUPERVISED LEARNING FOR SPLITTING COUNTRIES IN THE WORLD IN TERMS OF POSITION AGAINST SPREAD EFFECT

Alejandro Vaca Serrano

15 de enero de 2019

1. INTRODUCTION

2. Visualización Multivariante.

Transformaciones.

3. K-MEANS CLUSTERING

Using Euclidean Distance

Using Mahalanobis Distance

4. PAM

Using Euclidean Distance

Using Manhattan Distance

Using Mahalanobis Distance

5. Model-Based Clustering (Gaussian Mixtures)

6. Hierarchical Clustering.

Hierarchical Clustering with Euclidean Distance

Hierarchical Clustering Using Manhattan Distance

Hierarchical Clustering Using Maximum Distance.

Hierarchical Clustering Using Minkowski Distance

7(bis). Advanced Hierarchical Clustering.

Using Euclidean Distance.

Using Manhattan Distance.

Using Canberra Distance.

Using the “Complete” Method instead of Ward.D2.

7. Comparing The Best Clustering Models.

Comparación K-Means vs PAM

Comparación K-Means vs Hierarchical Clustering

Comparación PAM vs Hierarchical Clustering

Comparación Model-Based Clustering vs. Hierarchical Clustering.

COmparación K-Means vs. Model-Based Clustering.

Comparación PAM vs. Model-Based Clustering

8. Conclusiones